function [sigma_t, rho] = densityFunc(S, T, P)
  % Function for calculating density ala SBE Derive
  % S = salinity in PSU, T = temperature in deg C ITPS-68, P = pressure in bar
  %
  % Author: David Goldsmith, Wash. State Dept. of Ecology, dgol461@ecy.wa.gov
  % Release date: N/A

  % Constants from SBE Data Processing User Manual Revision 7.21d, 4/19/11,
  % Appendix V: Derived Parameter Formulas.
  % Implementation verified 4/12/2012
    B = [8.24493e-1; -4.0899e-3; 7.6438e-5; -8.2467e-7; 5.3875e-9];
    C = [-5.72466e-3; 1.0227e-4; -1.6546e-6];
    D = 4.8314e-4;
    A = [999.842594; 6.793952e-2; -9.095290e-3; 1.001685e-4; -1.120083e-6; 6.536332e-9];
    F = [54.6746; -0.603459; 1.09987e-2; -6.1670e-5];
    G = [7.944e-2; 1.6483e-2; -5.3009e-4];
    I = [2.2838e-3; -1.0981e-5; -1.6078e-6];
    J = 1.91075e-4;
    M = [-9.9348e-7; 2.0816e-8; 9.1697e-10];
    E = [19652.21; 148.4206; -2.327105; 1.360477e-2; -5.155288e-5];
    H = [3.239908; 1.43713e-3; 1.16092e-4; -5.77905e-7];
    K = [8.50935e-5; -6.12293e-6; 5.2787e-8];
  % Begin calculation
    T2 = T.*T;
    T3 = T.*T2;
    T4 = T.*T3;
    T5 = T.*T4;
    T = [ones(size(T)) T T2 T3 T4 T5];
    S(S <= 0) = 1e-6;
    S32 = S.^1.5;
    sigma = T * A;
    sigma = sigma + (T(:,1:5) * B) .* S;
    sigma = sigma + (T(:,1:3) * C) .* S32 + D*S.*S;
    sigma_t = sigma - 1000;
    kw = T(:,1:5) * E;
    aw = T(:,1:4) * H;
    bw = T(:,1:3) * K;
    k = kw + (T(:,1:4) * F) .* S + (T(:,1:3) * G) .* S32;
    k = k + (aw + (T(:,1:3) * I) .* S + J * S32) .* P;
    k = k + (bw + (T(:,1:3) * M) .* S) .* (P .* P);
    val = 1 - (P ./ k);
    sigma(val~=0) = sigma(val~=0) ./ val(val~=0);
    rho = sigma; 
end